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In his seminal work in the 1970s Robert May suggested that there was an upper limit to the 
number of species that could be sustained in stable equilibrium by an ecosystem. This deduction 
was at odds with both intuition and the observed complexity of many natural ecosystems. The so- 
called stability-diversity debate ensued, and the discussion about the factors making an ecosystem 
stable or unstable continues to this day. We show in this work that dispersal can be a destabilising 
influence. To do this, we combine ideas from Alan Turing’s work on pattern formation with May’s 
random-matrix approach. We demonstrate how a stable equilibrium in a complex ecosystem with 
two trophic levels can become unstable with the introduction of dispersal in space. Conversely, we 
show that Turing instabilities can occur more easily in complex ecosystems with many species than 
in the case of only a few species. Our work shows that adding more details to the model of May 
gives rise to more ways in which an equilibrium can become unstable. Making May’s simple model 
more realistic is therefore unlikely to remove the upper bound on complexity. 


It was once a widely-held belief in ecology that greater 
ecosystem complexity promoted stability [1-6]. A priori, 
this is an intuitive proposition. Surely the more species 
and connections there are in a food web, the less sensitive 
the network ought to be to perturbation. 

Providing a firm counterpoint to this view, Robert May 
used a simple statistical model [7, 8] to argue that in- 
creasing the number of species in an ecosystem could in 
fact reduce stability. By analysing the eigenvalues of a 
randomly constructed community matrix, May deduced 
the following criterion for stability [7] 


o NC <1. (1) 


In this criterion g? is the variance in the inter-species 
interactions, N is the number of species and C is the con- 
nectance (the probability that any given pair of species 
interact with one another). May’s result shows that sta- 
bility is decided by the product c = o? NC. We will call 
this quantity the ‘complexity’ of the ecosystem. 

May’s model suggests that more complex ecosystems 
can become unstable. For a fixed variance of interactions, 
there is an upper bound to the number of species and 
food web connections that the ecosystem can sustain. 
This idea quickly became controversial, and May’s work 
sparked the so-called complexity-stability (or diversity- 
stability) debate, which continues to this day [5, 6, 9]. 

Since the 1970s, the discussion around stability has 
been made more precise and subtle. It is now understood 
that there are a number of senses in which an ecosystem 
can be unstable, and indeed there are a number of ways 
one can define diversity [5, 10]. An ecosystem can be 
unstable with respect to, for example, the introduction 
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of new species, the extinction of existing species, or envi- 
ronmental changes. May’s work addresses stability with 
respect to fluctuations in species abundance. 

In order to understand the influence of the many as- 
pects of real ecosystems on stability, May’s fairly austere 
model has since been augmented and improved upon. 
Features not captured by May’s initial model include 
food-web structure (e.g., trophic levels, modularity, and 
nestedness) [8, 11-13], the feasibility of the equilibrium 
[14, 15], variation in interaction strength [16, 17], and 
variability of the environment and of species susceptibil- 
ity to environmental change [18-21]. In many models 
of complex ecosystems only the total abundance of each 
species is considered, without appreciating how the mem- 
bers of that species are distributed in space [7, 8, 11-21]. 
In such models there is no notion of space, and hence no 
dispersal. In this work we explicitly include the effects 
of diffusive dispersal in space, and study the effects on 
stability. 

While dispersal may intuitively be expected to be a 
homogenising and stabilising influence, Turing showed in 
his seminal work [22] that it can in fact destabilise a 
dynamical system. Such an instability has been stud- 
ied in meta-population predator-prey models with small 
numbers of species [23-26]. We combine Turing’s idea 
with May’s random-matrix approach to show that a sim- 
ilar destabilising effect can be seen in models of complex 
ecosystems. 

In order not to obscure the key effects at work, we opt 
to modify May’s paradigmatic model sparingly. This al- 
lows us to highlight the consequences of the inclusion of 
dispersal. We suppose that the abundances of the species 
rest in a steady, homogeneous equilibrium. In order to 
study stability, we examine the Jacobian matrix govern- 
ing perturbations about this equilibrium. Like May, we 
ask what statistical properties are required of the Ja- 
cobian matrix in order for the ecosystem to return to 
equilibrium when perturbed. 

Unlike May however, we allow for different trophic 
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The structure of the stability matrix M,. The matrix is composed of three parts: a diffusion matrix D, a self- 


interaction matrix d, and an interaction matrix A with entries from a Gaussian distribution. Each matrix is split into blocks 
due to the separation into two trophic levels — we use the subscript u to denote prey species and v for predator species. The 
approach can be generalised for greater numbers of trophic levels. The figure illustrates the content and structure of the blocks. 
The stability of the non-spatial ecosystem is described by the matrix for q = 0, or equivalently by setting Du = Dy = 0 (see 


Methods, and Section S1 in the Supplementary Information). 


levels, with statistical differences between them. It is 
the combination of dispersal and trophic structure which 
gives rise to the Turing-type instability. For the sake 
of mathematical simplicity, we confine our model to 
only two broad trophic classifications: predator and prey 
species. Our approach can be generalised to more com- 
plicated food web structures. 

We now postulate the form of the Jacobian matrix cen- 
tral to our problem, M,. The elements of this matrix de- 
scribe how spatial disturbances of wavelength A = 27 /q 
in the abundance of one species affect the abundances of 
the other species; q is known as the wavenumber [27]. 

Because of the separation of the population into 
trophic levels, the matrix M, has a block structure where 
each block has different statistics. The matrix is com- 
prised of three terms: a diffusion term, an intra-species 
interaction term and an inter-species interaction term. 
We write 


M; =- D-d+A. (2) 


The diffusion coefficients for prey and predator species 
are Du and D, respectively (Fig. 1). The interaction 
matrix A is modelled as having elements drawn from a 
Gaussian ensemble. Further details on the structure of 
M, and how one arrives at this form are given in Fig. 1 
and in the Methods section. 

The problem of analysing stability reduces to finding 
the eigenvalue spectrum of the matrix M4. If the eigen- 
values of this matrix all have negative real parts, then 
the equilibrium is stable with respect to disturbances of 
wavenumber q. Else, it is unstable. In order for the equi- 
librium to be stable on the whole, therefore, all eigenval- 
ues of M; must have negative real parts for all values of 


q. 
If the number of species in the ecosystem is large, then 
the eigenvalue spectrum is dependent only on the statis- 
tics of M, and not on its specific entries. Using random- 
matrix theory and ideas from statistical physics, we are 
able to deduce a mathematical expression for the support 
of the eigenvalue spectrum of M,. That is, we can find 
the region in the complex plane in which the eigenval- 
ues sit and, most importantly, whether or not they have 
positive real parts. Examples are shown in Fig. 2. 

With this analytical approach we can calculate what 
properties of M, make the equilibrium unstable. Thus, 
we can deduce how May’s upper bound on ecosystem 
complexity is modified by the inclusion of dispersal and 
trophic levels. Most crucially, we show that equilibria 
which would be stable without spatial effects can be 
destabilised by dispersal. We find that this dispersal- 
induced instability is possible not only in a linear model, 
but also in a non-linear system where the equilibrium is 
arrived at dynamically, and hence feasible by construc- 
tion. 


Results 

Eigenvalue spectra. We show some example eigen- 
value spectra of the matrix M, in Fig. 2. The vast ma- 
jority of the eigenvalues group into a ‘bulk’ region, with 
the exception of a few outliers. These outliers cannot 
be ignored — the excursion of even one eigenvalue across 
the imaginary axis to the positive real side makes the 
equilibrium unstable. 

Using the statistical properties of M, we are able to 
calculate mathematically the bulk regions to which most 
of the eigenvalues are confined, and the locations of any 
outliers. In Figs. 2 (a)-(c) we show that these calcu- 
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FIG. 2. Eigenvalue spectra of the stability matrices in Fig. 1. In May’s original model, the eigenvalues all lay uniformly within 
a circle in the complex plane. In our model the circle is warped and split into more complicated shapes. A small number of 
outlier eigenvalues now stray from the bulk. Eigenvalues of computer-generated random matrices M4 are shown as blue crosses, 
and compared to theoretical predictions for the boundary of the bulk region (red solid line). Predictions for the outliers are 
show as open red circles. For q = 0 all of the eigenvalues have negative real part (panel (a)). The equilibrium is stable in the 
non-spatial ecosystem. For disturbances with larger wavenumber q, one of the eigenvalues crosses the real axis (panel (b)). The 
equilibrium is unstable against such perturbations. If the wavenumber q is larger still the most unstable eigenvalue returns to 
the negative half-plane (panel (c)). This is characteristic of a Turing instability [27] — the equilibrium is unstable with respect 
to disturbances of a finite range of wavelengths. This is shown in panel (d), where we plot the real part of the most unstable 
eigenvalue as a function of q. The equilibrium is unstable against perturbations of wavenumber q whenever Re[wmax] > 0. The 


green triangles mark the wavenumbers from panels (a)—(c). 


lations agree very well with the spectra of computer- 
generated random matrices. We can therefore predict 
what community properties lead to stability or instabil- 
ity. 

As Fig. 2 demonstrates, it is possible to find cir- 
cumstances under which the model community is desta- 
bilised by the inclusion of dispersal. Fig. 2 (a) shows 
the eigenvalue spectrum for the model without spatial 
effects (q = 0, or Du = Dy = 0, see Methods). All 
eigenvalues in panel (a) have negative real parts and we 
conclude that the equilibrium is stable for the non-spatial 
model ecosystem. In Fig. 2 (b) we take into account dis- 
persal and show the eigenvalue spectrum for a non-zero 
wavenumber. All other parameters are the same as in 
panel (a). An outlier eigenvalue strays over the imagi- 
nary axis, demonstrating that the equilibrium is now un- 
stable. For perturbations with higher wavenumbers the 
outlier then returns to the negative half-plane (panel (c)) 
— the equilibrium is stable with respect to perturbations 
of higher wavenumber. 


Fig. 2 (d) shows the real part of the most unstable 
eigenvalue as a function of the wavenumber q. We see 
that the equilibrium is unstable against perturbations in 
a finite band of wavenumbers. In Turing’s original work 


on chemical reaction systems [22], this signalled the for- 
mation of stable periodic patterns. The exact shape of 
these patterns is usually determined by non-linearities in 
the differential equations describing the reactions. Our 
model is valid only in the vicinity of the supposed equi- 
librium, and similar to May [7] we have not specified the 
nature of any non-linearities. We therefore do not spec- 
ulate for now about what might happen after the system 
has departed from the fixed point. We merely point out 
here the dispersal-induced instability of the equilibrium 
about which we have linearised. 


Modifiying May’s bound: stability with and with- 
out dispersal. In order to further appreciate the ef- 
fect that the inclusion of dispersal has on stability, we 
first consider the conditions under which the non-spatial 
ecosystem becomes unstable (see Methods). This en- 
ables us to study how May’s bound on the complexity 
c changes for our model, which has distinct predator and 
prey species. A stability plot is shown in Fig. 3 (a). 
The horizontal axis shows the average degree of predation 
p = CN hvu (see Methods), the vertical axis is the com- 
plexity parameter c. The solid line indicates the upper 
bound on the complexity: below the line the equilibrium 
is stable, above this line it is unstable. 
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FIG. 3. Dispersal as a destabilising influence. Panel (a) shows the stability diagram for the non-spatial ecosystem model with 
two trophic levels. The solid line is the upper bound on complexity c for the equilibrium to be stable. Increasing complexity 
leads to instability. There is a lower bound on the amount of predation p required for the community to be stable (if p is 
too small, then the equilibrium is unstable even for small values of the complexity). Stability for the spatial ecosystem with 


dispersal is shown in panel (b). 


In the blue region the equilibrium is stable. In the yellow region the equilibrium would be 


stable in a non-spatial model, but dispersal in space induces instability. Crossing from the blue into the yellow region in panel 
(b) the ecosystem undergoes a Turing instability. In the red region in (b) the equilibrium is unstable both in the non-spatial 


and in the spatial ecosystem. 


We see from Fig. 3 that greater predation p increases 
the amount of complexity c that can be sustained in sta- 
ble equilibrium by the ecosystem. Notably, in order to 
have stability at all there is a lower bound on the preda- 
tion parameter p. 

If we now include dispersal, the stability diagram 
changes (Fig. 3 (b)). In particular the upper bound 
on the complexity can become lower than in the non- 
spatial system. This is because a new type of instability 
is now possible — the Turing instability. Thus, there are 
instances in which the model is stable without dispersal, 
but unstable when dispersal is introduced (yellow area in 
Fig. 3 (b) labelled ‘dispersal-induced instability’). 

There are no situations in which an unstable equilib- 
rium is stabilised by dispersal (see Section S5 in the Sup- 
plementary Information). Dispersal is a purely destabil- 
ising influence in our model. 


How does complexity affect the Turing instabil- 
ity? So far we have concerned ourselves with the effects 
of dispersal on complex ecosystems. We now ask the re- 
verse question: Spatial instability and pattern formation 
have been found in ‘simple’ models of ecosystems with a 
small number of species [28-31]. What are the effects of 
complexity on this Turing instability? 

Turing instabilities in simple systems typically occur 
when the diffusion coefficients of the activator and in- 
hibitor components are quite disparate — this is known 
as the ‘fine-tuning issue’ with the Turing mechanism 
[32, 33]. The activating components in our system are 
the prey species, predators play the role of the inhibitors. 


We would like to know what the effects of complexity are 
on the threshold ratio D,/D,, at which the Turing insta- 
bility occurs. To answer this question, we compute this 
threshold for different values of the complexity parameter 
c. 

The case c = 0 can be achieved by setting the width ø 
of the Gaussian distribution for the elements of the ma- 
trix A to zero. All matrix elements within each block are 
then identical to each other. This means that all preda- 
tor species become identical to one another, and similarly 
there is no distinction between different prey species. The 
model reduces to a simple two-species predator-prey sys- 
tem. Increasing c from zero introduces heterogeneity be- 
tween the species. 

We find that complexity can lower the ratio of diffusion 
coefficients required for Turing instability (Fig. 4). For 
more complex ecosystems the Turing instability therefore 
sets in more easily. So, not only can complexity decrease 
the stability of non-spatial model ecosystems, it can also 
reduce stability in spatial models. Conversely, increas- 
ing the ratio of diffusion coefficients D,/D, reduces the 
complexity that can be sustained in stable equilibrium. 
As can be seen in Fig. 4, the upper bound on c is lower 
when disparity of dispersal rates is large. 


Spatial instability in a non-linear model with 
complexity. The linear analysis we have focused on so 
far, although informative, has drawbacks. It only deals 
with the dynamics in the vicinity of a homogeneous equi- 
librium, and it tells us nothing about how the ecosystem 
behaves in the long-run if the equilibrium is unstable. 
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FIG. 4. Complexity reduces the ratio of diffusion coefficients required for Turing instability. In the yellow region the community 
is unstable due to Turing instability. A high value for the ratio of diffusion coefficients D, /Du reduces the complexity that 
can be sustained in stable equilibrium. For sufficiently high complexity the non-spatial community becomes unstable (red area 
to the right of the vertical line). The equilibrium of the spatial ecosystem is then also unstable, irrespective of the diffusion 


coefficients. 


Further, one could object that the linear model is some- 
what contrived and that it does not capture how a ‘real’ 
ecosystem constructs itself in equilibrium. 

We now present simulation data of a complex ecosys- 
tem obeying non-linear Levin-Segel-type dynamics [26] 
(see Methods). In Fig. 5 we demonstrate that a 
dispersal-induced instability can occur in this model as 
well. Panel (a) shows a realisation of the dynamics with- 
out dispersal. The model ecosystem converges to an equi- 
librum composed only of surviving species. By definition 
this is a feasible equilibrium [34]. 

Fig. 5 (b) shows the same model (with the same in- 
teraction matrix A), but now with dispersal. The abun- 
dances do not settle in the long run. Instead, they display 
quite erratic behaviour. We stress that this is different 
to the typical behaviour seen in two-species systems with 
a Turing instability. In such simple systems, the abun- 
dance of each species typically settles down eventually. 
The fixed-point value varies across space, generating a 
periodic pattern. The complex nature of the interactions 
in Fig. 5 (b) leads to more complicated dynamics known 
as diffusion-induced chaos [35, 36]. 

That being said, at any point in time one can take a 
snapshot of the spatial profile of the species abundances. 
An example is shown in Fig. 6. One finds some spatial 
structure for the abundances of prey species (green lines), 
but it differs from classic Turing patterns, which are typi- 
cally more periodic and regular. We note that the quickly 
diffusing predator species (red lines) have more smoothly 
undulating spatial profiles than the slowly diffusing prey 
species in this example. 


Discussion 
The random-matrix approach to modelling ecosystems 


has developed substantially since May’s original work. 
We contribute to this development with a model of a 
complex ecosystem which includes different trophic levels 
and dispersal in space. These additional features change 
May’s bound on complexity. 

How the bound changes tells us about the influence of 
the new model components on stability. For example, we 
find that predation acts as a stabilising influence (Fig. 3, 
and Section S5 of the Supplementary Information). 

Inspired by Turing’s mechanism for pattern formation 
[22], we show how dispersal can be a destabilising factor 
in a complex ecosystem. An equilibrium which would 
be stable in a non-spatial model can become unstable 
through dispersal. Conversely, we observe that increased 
complexity can lower the threshold for the diffusion co- 
efficients required for a Turing instability. So, not only 
does complexity reduce stability in a non-spatial model 
as was May’s conclusion [7], it also destabilises spatial 
models with dispersal. 

This is interesting especially in light of the so-called 
‘fine-tuning’ issue with the Turing mechanism. The ratio 
of diffusion coefficients required for a Turing instability is 
usually large, making it hard to find experimental exam- 
ples of Turing pattern formation [27, 33, 37, 38]. Based 
on our findings, we speculate that Turing patterns may 
be easier to observe in complex dynamical systems. 

To demonstrate that the Turing instability can also be 
seen in model systems where the equilibrium is arrived at 
more naturally, we performed simulations of a complex 
ecosystem with Levin-Segel dynamics (Fig. 5). We found 
that this model also has a dispersal-induced instability. 
Because of the complex nature of the interactions the 
system does not reach a stable patterned state, normally 
characteristic of Turing instabilities in models with fewer 
species. Instead, one sees persistently volatile dynamics 
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FIG. 5. Dispersal can induce instability of feasible equilibria. The prey species abundances in the non-spatial Levin-Segel 


ecosystem in panel (a) converge to a stable feasible equilibrium. 


In panel (b) we simulate the ecosystem with the same 


paramaters as in (a), but allowing for dispersal in space. Volatile behaviour is found in the spatial ecosystem, even though the 


non-spatial system would approach a stable equilibrium. 


-— 


1" 


HAN 


JAN 


FIG. 6. Distribution of species in space as a result of dispersal-induced instability. The figure shows a snapshot of all N = 100 
species abundances in a complex Levin-Segel ecosystem. The model parameters are the same as in Fig. 5 (b), and given in full 
in the Supplementary Information. Predator species are shown as red lines, prey species as green lines. The most abundant 
species varies from place to place. This is structure is not a static pattern as in conventional reaction-diffusion systems with 


Turing instability. Instead the abundances change with time. 


known as diffusion-induced chaos [35, 36]. 


Our notion of dispersal as a destabilising influence is 
in contrast to a recent study by Gravel et al. [39] who 
assert that dispersal can in fact stabilise a complex meta- 
ecosystem. One key difference between this work and 
ours is the way in which dispersal is implemented. In Ref. 
[39], the ecosystem is organised into patches, and disper- 
sal occurs from any patch to any other patch with equal 
probability. In other words there is all-to-all connection 
between the patches, and therefore no sense of spatial 
arrangement of these patches. In our model, dispersal 
is local and mediated via a standard diffusion term. We 


regard local dispersal as a more realistic way to model 
species movement. The fact that our findings are differ- 
ent from those in [39] indicates that the way one models 
dispersal can have quite an impact on its effects on sta- 
bility. We note that the theory we have developed can 
also be used to predict stability in meta-ecosystems with 
local dispersal between neighbouring patches (see Section 
S1 of the Supplementary Information). 


One criticism levelled at May’s model is that it is too 
simple and that perhaps through the inclusion of further 
aspects of natural ecosystems, the upper bound on the 
complexity would be alleviated. Our results do not sup- 


port this hypothesis. Like May, we also find that there is 
always an upper limit on the complexity c = o? NC that 
an ecosystem can stably sustain. Other recent studies us- 
ing random-matrix approaches arrive at similar conclu- 
sions. For example, Allesina and Tang take into account 
more realistic food-web structures and still find upper 
limits on the number of interconnected species [12, 13]. 
May’s result therefore generalises to models capturing 
more aspects of ‘real’ communities in ecology. This pre- 
diction is supported by the observation of ‘diversity reg- 
ulation’ in some ecosystems [40—44]. 

A final observation that we wish to convey is that mak- 
ing models for complex ecosystems more detailed intro- 
duces the opportunity for new types of instability. In 
May’s original model, for example, the mean of the com- 
munity matrix elements was zero. As a consequence, any 
one species is equally like to benefit or suffer from the 


presence of another species. Mathematically, all eigenval- 
ues then reside within one bulk region, and it is this bulk 
region that determines stability. Mutualism can be intro- 
duced through interaction coefficients which are positive 
on average, and competition through a negative average 
interaction [12]. This leads to additional outlier eigenval- 
ues, which can make an equilibrium unstable even though 
it would otherwise be stable. Introducing trophic levels 
can generate complex-conjugate pairs of outliers (Fig. 2 
a)), allowing for more opportunity for instability (Sec- 
tion S5 of the Supplementary Information). Dispersal, 
finally, leads to the possibility of a Turing-type instabil- 
ity. Overall, adding more details to the model of May 
tends to give rise to more ways in which an equilibrium 
can become unstable. Making May’s simple model more 
realistic is therefore unlikely to remove the upper bound 
on complexity. 


— Methods — 


Linear model. We imagine that we find the ecosystem at an homogenenous equilibrium. Our model is concerned 
with the dynamics of small perturbations of the species abundances about this fixed point. The stability of the 
homogeneous fixed point is determined by whether or not these perturbations decay or increase with time. We write 
u;(x,t) and v,;(«,t) for the perturbations of the prey and predator species abundances respectively at position xz and 
time t. These are the deviations away from the fixed point. There are N, prey species and N, predator species with 
N = N, + N, species in total. We define the constants yu = N,/N and y = N,/N. 

We assume that prey species diffuse at rate D, and predator species at rate D, in a spatially homogeneous 
environment. Similar to May [7], we also imagine that all species have negative self-interaction. This is a necessary 
condition for stability and reflects intra-species competition. The probability that a particular pair of species interact 
with one another is C. This parameter is known as the ‘connectance’ [7]. The effect of a change in the abundance of 
species j on species 7, where j belongs to trophic level 8 and i belongs to a, is AGP la, 8 € {u,v}]. The linearised 
reaction-diffusion equations thus take the following form 
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If species i and j are non-interacting, ASP = a = 0. If the two species interact, then A and Ave are drawn 


from a joint Gaussian distribution with means Age = Uap and Ate = Ugo. The elements in the random matrix have 


variance g? 


(ASP — pag)? = 0, (4) 


and they are correlated according to 


(ARP = Hap) AS? = pogi) = Tagapo’. (5) 


All other correlations are set to zero. Each of the model parameters can be interpreted ecologically. We define 
and interpret the complexity c = CNo? in the main text. The interaction mean Huu indicates the degree to which 
different prey species cooperate (if puu > 0) or compete (if Huu < 0). The coefficient jz», has a similar interpretation 
for predators. The means of the off-diagonal blocks fu» < 0 and Hyu > 0 indicate the degree to which (on average) 
prey species suffer and predator species gain from predator-prey interactions. The parameters Pygq’g describe the 
correlations between interaction coefficients. The only non-zero entries are taken to be Ty = Tyuuu, Lo = Povyy and 


Tuv = Tuvvu- That is, only elements which are diagonally opposite one another in A are correlated. A positive value 
of Taga'g' indicates that if one species benefits more than average from an interaction, the other species involved does 
so as well. The opposite is true if [ggq7g’ is negative. 

Taking the Fourier transform with respect to the spatial coordinate x of Eq. (3), we arrive at dynamical equations 
(see Supplementary Information Section S1) for disturbances of wavenumber q in the abundances of the various 
species (the wavenumber is related by q = 27/X to the wavelength A). We denote the combined vector of the Fourier 
transforms of species abundances by X, = (--- ,u(q,t),--- ,0;(¢,t),--+), and arrive at the more compact matrix 
equation 


X, =M,Xq. (6) 


The matrix entry (M,)¢" tells us what the effect of a disturbance of wavelength q in species j (belonging to trophic 
level 3) is on species i (belonging to trophic level a). 

The vector X, has dimension N = N, + Ny and is arranged such that the first N, elements are the Fourier- 
transformed abundances &;(q, t) and the last Ny elements are the 0;(q,t). The matrix M, is depicted in Fig. 1. It has 
a block structure due to the separation into trophic levels. Its three contributions are a diagonal diffusion matrix, a 
diagonal self-interaction matrix and a random interaction matrix, whose variance and correlations are given in Eqs. (4) 
and (S2). The indices a and 8 correspond to the different blocks and i and j correspond to the position within the 
block. If the matrix M, has eigenvalues with positive real parts for any value of q > 0, the disturbances u; and v; 
will grow with time, indicating an unstable equilibrium. 

By setting q = 0 in Eq. (6), one recovers Eqs. (3) with the diffusion term removed. Focusing on q = 0 in our model 
(or equivalently setting D,, = D, = 0) thus allows one to study the stability of a non-spatial model ecosystem with 
two trophic levels. See also Sections S1 and $5 in the Supplementary Information. 

We note that the same form for the matrix M, can be found for a model of a meta-ecosystem (with dispersal 
between nearby patches instead of in continuous space), see Section S1 of the Supplementary Information. 

The values of the model parameters used in the figures are given in full in Section S6 of the Supplementary 
Information. 


Calculation of the boundary surrounding the bulk of the eigenvalues. The vast majority of the eigenvalues 
of the random matrices M, reside within one or two ‘bulk’ regions of the complex plane. To determine stability we 
need to know if there are bulk eigenvalues with positive real parts. Identifying the boundaries of the regions containing 
the eigenvalues is sufficient for this purpose. Our calculation uses methods originally developed in condensed matter 
physics, and follows lines similar to those of [45, 46]. 

We write w = wz, + iwy for the eigenvalues of M,. The general expression for the boundary surrounding the 
bulk of the eigenvalues (wy as a function of wg) is given by the simultaneous solution of the following equations (see 
Supplementary Information Eq. (S81) ) 


bac Xal — : = 0, 


— (wz + iwy + da + @Da)Xa + cX Tap yexaXs +1=0. (7) 
B 


We note the free index a in the second of these equations (œ € {u, v}). This is therefore a system of three coupled 
equations. One first eliminates the auxiliary variables Xa, and then expresses wy in terms of w,. This results in the 
red curves in Fig. 2. 

The solution simplifies in several special cases, which we exploit to provide explicit stability criteria analogous to 
May’s bound (Section S5 of the Supplementary Information). 


Calculation of outliers. In addition to the bulk eigenvalues the stability matrix can have isolated outlier eigenvalues. 
If any of these outliers have positive real part, the equilibrium is unstable. Their position in the complex plane is 
calculated along the lines of [47]. Details can be found in the Supplementary Information. The outlier eigenvalues 
are given by the complex values w satisfying the following equation (see Supplementary Information Eq. (S82)) 


1 
YuNC' Huu — | [oN Cu a Xo (w) = (NC)? Yu Yvhuvhvu =0, (8) 


Xv(w 


= 
Xu(w) 


The auxiliary quantities Xa(w) in this relation satisfy 


—1 = —(w + du + q Du)Xu =F Vuyuxe + uv YuXoXus 


1 = —(w + dy +° Do)Xo + Cur yuxuxe + Vow; (9) 
subject to the condition 
2 1 
D elxa(w)? < =. (10) 


Eqs. (S20) and (S21) need to be solved simultaneously, subject to Eq. (S22). If there are no solutions then there 
are no outliers. In special cases the above expressions can be simplified, and explicit stability criteria can be found 
(Section S5 in the Supplementary Information). 


Finding the threshold for instability. Instability can occur in one of several different ways: (1) The bulk region 
of eigenvalues for q = 0 can cross into the positive half-plane; (2) One of the outlier eigenvalues for q = 0 can cross 
the imaginary axis; (3) An outlier eigenvalue for q Æ 0 can stray into the positive half-plane. We have not observed 
any circumstances under which the bulk crosses the imaginary axis for non-zero q where it does not for q = 0. The 
eigenvalues must have negative real parts for all q in order for the equilibrium to be stable in the spatial system. This 
includes g = 0. Cases (1) and (2) therefore indicate instabilities occurring both in the non-spatial and the spatial 
ecosystem. In case (3) the spatial system is unstable, but the non-spatial system remains stable. In each of these 
cases, the threshold for instability is found by identifying sets of parameters for which either the boundary for the 
bulk eigenvalues touches the imaginary axis (case 1) or where the outlier eigenvalues touch the imaginary axis (cases 
2 and 3). 

This can be done using our analytical results for the spectrum of eigenvalues (Section S5 of the Supplementary 
Information), leading to the results shown in Figs. 3 and 4. 


Simulating the non-linear model. Results in Figs. 5 and 6 are from a numerical integration of the Levin-Segel-type 
model [26], 


Ou; Oru; uu uv 
ae = Duga ti a— uit) Aur + > Ato; ; 
= keu jeu 
Ov; O70; , 
y =D aut Aant DD A (11) 
icx kev 


where a > Oisaconstant. To integrate these equations numerically, the diffusion terms are discretised. The integration 
is then carried out using the Runge-Kutta (RK4) method [48]. 
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— Supplementary Information — 


S1. MODEL CONSTRUCTION 
A. General setup: trophic levels, species interactions and diffusive dispersal 


We suppose that the total population consists of a sub-population of predator species and a sub-population of prey 
species, which are distributed in space. We assume that there are Nu = y_N different prey species and N, = YN 
predator species, where Yu + Y% = 1. N is the total number of species, N = N, + Ny 


We imagine that the system reaches an equilibrium and that the equations governing the dynamics of the species 
abundances, linearised about this equilibrium, take the following form 


Ou; po i 
p” “De a = a+ Ye At + Ye Ay vj, 
Ov; p? vj f 
Fe = DT — deny YAU +Y An s) 
The u;(x,t), i = 1,..., Nu, describe perturbations of the abundances of the prey species about the equilibrium at 
position x and time t, and similarly v;(a,t), j = 1,..., Ny, are perturbations of the predator abundances. 


We presume that the interactions within a species are competitive (such that d, > 0 and d, > 0), and that 
each species disperses via diffusion. The corresponding dispersal coefficients are D,, for prey, and D, for predators. 
Crucially, the prey and predator species can have different dispersal coefficients Du # Dy. In the model there are 
inter-population and intra-population interactions. These are described by the interaction coefficients Ag? , where 
a and § can takes values a, 8 € {u,v}, where u stands for prey and v for predator. For example, Aj” describes 
interactions between different prey species, and Ai? describes the effect that predator species j has on prey species 7. 


The interaction coefficients AY are random ean drawn from a Gaussian distribution with 


ap 
Ai; = Haß, 
(AS? — pap)? = 0°, 
(ASP — pap (AS — Harp) = Tagpap o’. (52) 


All other correlations are zero. We assume that the only non-zero Tages are Puuuu = Fu, Cuvvu = Pus and 
Povvv = Ty. With a probability 1 — C, species i and j are chosen not to interact with one another and thus Age and 


AG are set to zero. The parameter C therefore describes the ‘connectance’ (the probability that any given pair of 
species interact with one another). 

Each of the model parameters can be interpreted ecologically. We define and interpret the complexity c = CNo? in 
the main text. The interaction mean Huu indicates the degree to which two different prey species cooperate on average 
(if Huu is positive) or compete (if Huu is negative). The coefficient jy, has a similar interpretation for predators. The 
means of the off-diagonal blocks Huv < 0 and fy, > 0 indicate the degree to which (on average) prey species suffer 
and predator species gain from predator-prey interactions. The coefficients I'4gq’g: indicate statistical correlations 
between interaction coefficients. A positive value of [agq’g’ indicates a statistical tendency for one species to benefit 
more (than average) from an interaction if the other species involved does so as well. The opposite is true if Paga’g: is 
negative. The only non-zero entries are taken to be Ty, = Puuuu, Po = Vvvvy and Puy = Tuvvu. That is, only elements 
which are diagonally opposite one another in A are correlated [49, 50]. 


B. Fourier transform and stability matrix 


The perturbations u;(x, t) and v,;(z,t) in Eq. (S1) are functions of position (and time). Taking the Fourier transform 
with respect to the spatial coordinate x, we write 


: dtu (a 
= sm [are a( ); (S3) 
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and similar for ù; (q). The ù; (q), 0; (q) describe perturbations of wavenumber q, or equivalently, of wavelength A = 27 /q. 
Using this in Eqs. (S1) one obtains 


Nu Ny 
ùi = —( Du + duis +X ARI, + SO AMO, 
j=l 


k=1 
Nu Ny 

bj = —(°Dy + dy); +) Apts + DART, oe 
i=1 k=1 


We note that by setting q = 0 in Eqs. (S4) one recovers Eqs. (S1) with the diffusion term removed. From the definition 
of the Fourier transform Eq. (S3), we see these are then equations for the total u; and v; summed throughout space. 
This removes any sense of spatial distribution of species. Setting q = 0 therefore recovers the naive ‘non-spatial’ 
model where species dispersal is not taken into account. Eqs. (S4) can be written in a more compact matrix form 


X=M,X, (S5) 
where the vector X is of length N, and with 
M; = -7 D-d+A. (S6) 
The matrix M, is of dimension N x N. It is the sum of three terms: 


(i) The first term is —q?D, where the matrix D is of size N x N, diagonal and of the form 


a f Duln,, 0 
p= (24n 2), 6 


where Iyn, and 1y, are identity matrices of dimensions N, x Nu and N, x N, respectively. 
(ii) The second term is of a similar form 


_ { din, 0 
an(4 42), es 
and describes intraspecific interaction. 
(iii) The third term in M, is a random matrix with the Gaussian statistics given in Eq. (S2). We note that the 
matrix M, has a block structure; it is comprised of submatrices of dimensions N, x Ny, Nu X Ny, Ny X Ny and 
Ny x Ny. This is illustrated in Fig. 1 of the main text. 


The eigenvalues of the matrix M, tell us about the stability of the hypothetical fixed point about which we have 
linearised. More precisely, the fixed point is stable against perturbations of wavelength A = 27/q if and only if all 
eigenvalues of M, have negative real part. The fixed point is stable against perturbations of any wavelength only if 
this is the case for all q. That is to say, in order for the equilibrium to be stable, the matrix M, must not have an 
eigenvalue with positive real part for any q. 

We therefore wish to identify the eigenvalue spectrum of M,. We can then use this to deduce under what conditions 
the equilibrium is stable. The matrix My has N, + N, eigenvalues in total. The eigenvalue spectrum consists of a 
‘bulk’ of eigenvalues with some possible outliers (see Fig. 2 in the main text). The bulk contains the vast majority of 
eigenvalues. They are confined to a compact area in the complex plane. Both the boundary of the bulk region and 
the location of the outliers can be calculated based only on the statistical properties of the matrix M, (as opposed 
to the exact entries). Further, the boundary of the bulk region can be calculated independently of the outliers. We 
therefore split our calculation into two parts, and present the calculation of bulk eigenvalues in Section S2, and that 
for the outliers in Section S3. 


C. Correspondence with a model of a meta-ecosystem with patches 


We now briefly consider a model in which species hop between discrete patches, as opposed to diffusion in a 
continuous space. We label patches by an integer index £Z, and assume that patch £ is located at position Ax in space. 
That is to say, the spacing between two adjacent patches is Az. The continuous diffusion term in Eqs. (S1) is then 
replaced to yield 


z = 2D, XO [U Lyu (l Ar) — oU, Oui (lAr)] — duu (lAr) 


VA 
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53 uuu, (LAr) +e ip (Az), 


Ov;(CAr) _ A 1\n,.( pt _ ’ 0\u; x)| — d,v; r 
= 2D, X [06 l)u; (Az) — (l, &)v;(CAx)] — duv; (Az) 


7 
+S agra (€Ax) + 53 Azp (Az). (S9) 


(We have suppressed the dependence on time to keep the notation compact). The quantity (£, 0’) is a ‘hopping 
kernel’, describing the probability that an individual hops from patch ¢’ to patch £, given that a hop occurs. We 
suppose hopping is only between neighbouring patches, and occurs with equal probability. We therefore choose 


$ if €=£ +1, 
ot, 0) = l if €=f-1, (S10) 
0 otherwise. 


The constants D,, and D, in Eqs. (S9) describe the overall hopping rates (for prey and predators respectively) to any 
of the two neighbouring patches. If the migration between patches is diffusive (or just by dimensional arguments) 
these rates must be proportional to (Ar)~?, so we write D, = D,/(Ax)? and D, = D,/(Ax)?. 

Next, we take the discrete Fourier transform with respect to the spatial variable in Eqs. (S9) , 


q) = 5 ÀLA (LAT). (S11) 
4 

We obtain 

Dü a an 

a= 2D,,[6(q) — Iù; — d vie a thai tÀ at 

06; a 

a = 2D,[6(q) — 1]õ; — dt; + year i; + S.A, (S12) 
with 


plq) = cos (qAr). (S13) 


Provided the arrangement of patches is ‘dense’, i.e., that there are many patches per unit distance, we can assume 
Aa <1, and expand ¢(x) = cos (qAz) in powers of Ax. To lowest non-trivial order we find 


~ 1 
olz) x1- z147)’. (S14) 
Eqs. (S12) then become 
OU; 2 ~ i uu ~ Z 
ap = (P Du + du) XO Atttin + 55 AHD, 
k=1 j=l 

00; ~ vug, 
= = ~(q2D, + dajë; yay +S Ae (S15) 


and we thus arrive at Eqs. (S1). 

This shows that any stability criteria derived for the ecosystem with diffusion in continuous space also applies to a 
model with discrete patches and hopping between adjacent patches. The approximation is valid provided the distance 
between two patches in space is small (Az < 1), and we have assumed the standard scaling D, = D inf (any, Dy = 
D,/(Az)?. It is also possible to analyse the model with discrete patches, and to find the eigenvalues as a fancdod of 
q without performing this approximation. One would still find the possibility for Turing instability (as demonstrated 
for a different system in [33]). 
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S2. BOUNDARY OF THE BULK OF THE EIGENVALUE DISTRIBUTION 


In this section, we broadly follow the structure of similar calculations in [45, 46, 51]. 


A. Preliminaries — setting up the calculation 


For an ensemble of N x N random matrices m we consider the so-called resolvent, defined by 


G(w) = <r [wy —m)7!]. (S1) 


The overbar represents an average over the ensemble of matrices. The expression 1, denotes the identity matrix of 
size N x N. 
We note the following identities from general linear algebra (for an N x N matrix B): 


(i) Tr [P~! BP] = Tr [B], for any invertible N x N matrix P; 

(ii) If B is invertible, then the eigenvalues of B~~ are reciprocals of those of B; 

(iii) The trace is the sum of the eigenvalues; 

(iv) Adding the identity matrix to a matrix increases each eigenvalue of this matrix by 1; 


(v) Multiplying a matrix by a constant scales its eigenvalues by that same constant. 


Using these, one can show that 


GW) =7> (82) 


where the A are the eigenvalues of m. In the limit N — oo, this sum can be written as an integral over the complex 
plane 


2) PA) 
G (w) = fa aL, (S3) 
where p(A) is the density of eigenvalues. 


Bearing in mind that G(w) is a complex quantity, we perform a contour integral around a closed path C, which we 
presume contains none of the poles of G, 


aa fe G(w) 
2ri Je 


1 T T 
NX mi | 
2 1= [ dA pA), (S4) 


In this expression S denotes the region bounded by C, and we have used the residue theorem to evaluate the contour 
integral. 
We next note the complex version of Gauss’ law (letting w = x + iy) 


[ee Ea | - z [aw ($5) 


Defining the real-valued potential ®(w) via 


: l 


A E jimea A, (S6) 
Ox 


and using Eqs. (S4) and (S5), along with the fact that the region S is arbitrary, we obtain Poisson’s equation 


V°® = —4rp. (S7) 
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In analogy with electrostatics, the resolvent G(w) can be thought of as being related to the ‘electric field’ with 
components Ey = 2Re G and Ey = —2Im G [46]. Therefore we can obtain the distribution of eigenvalues, if we can 
find the potential ®. 

The potential ® can be related directly to the matrix m via 


(w) = -cm det [(1yw* — mT) (Lyw — m)]. (S8) 


This can be seen to agree with Eqs. (S6) and (S2) once one realises that det (Iw — m) = []\(w — A). 
We note the following about this potential and its relationship with the eigenvalue density. By definition, the 
resolvent (or average Green’s function) is given by 


o 1 1 _ (w, w*) 
G(x,y)= wit Fe = =| Iw (S9) 
Indeed, we also have 
O®(w,w*) 
x = 2 1 
G* (x,y) Ju (S10) 
The resolvent is related to the eigenvalue density by 
EG, e | xpi eae )) ope a (S11) 
Ox Oy Ow* 


The left-hand side of this equation is familiar from the Cauchy-Riemann equations for analytic functions. From this 
we deduce that if G(x, y) is an analytic function in a particular region (i.e. we can write it in terms of w only), the 
eigenvalue density must be zero in this region. 


B. Carrying out the average over the ensemble of random matrices 
Ensemble of random matrices 


Adding a low-rank perturbation to a large random matrix is known not to change the bulk of the eigenvalue 
spectrum [12, 47, 52]. This means that we can assume that the mean of the distribution from which the matrix 
elements are drawn is zero. The mean can have an effect on outlier eigenvalues, and we will therefore restore it for 
the calculation of the outliers in Section 53. 

For the time being we consider a Gaussian random matrix with statistics as follows 


i =0, 
2 
apja 7 
(mi; ) =W’ 
—; ap ~— Ta Bargr o? 
B Bi _ + aBa’B 
ni = Pate 
me = —da, (S12) 


The matrix m is similar to M, (the matrix we focus on in the main paper), but it is different in a number of ways: 

First, we note that we have re-scaled the matrix elements relative to the matrix M, in the main paper — the 
variance and co-variances are now proportional to 1/N. This approach is standard in the theory of random matrices, 
its objective is of a technical nature. The re-scaling allows us to carry out the calculation in the limit N — oo. 
In practice the results derived in this way are often accurate also for finite N. The scaling is undone by replacing 
a? — No? in the result. For a more detailed discussion, see Section $4. 

Further, we are assuming all-to-all interaction (C = 1) for now between the species. We discuss how choosing a 
general connectance 0 < C < 1 affects the result in Section S4. 

Finally, the wavenumber q is set to zero here. The result for non-zero q can obtained by making the replacement 
da — da + q?Dqa, see also Section S5 B. 


15 
Replica method and Hubbard-Stratononich transformation 


In order to calculate the spectrum of the ensemble of random matrices we need to compute 
Indet [(w*1_y — m7)(wiy —m)]. We recall that the overbar stands for the average over the randomness in m. In 
order to calculate averages of logarithms, one can use the so-called replica method [46, 52-54], based on the identity 

¥—! This allows one to calculate y” instead of the more intricate average In y. The object y” in turn 
represents an n-fold ‘replicated’ system in statistical physics [46, 52-54]. In the context of the present problem one 
finds that the n replicas decouple in the limit of large N, similar to what was observed in [45, 46]. This means that 


In det [(1yw* — mT) (Lyw — m)] = Indet [(1yw* — m™)(1yw — m)]. (S13) 


Using a Gaussian integral we can write the determinant in Eq. (S13) in the form 


det [(yw* —m™)(yw — m)] * = i: II (£) 


xexp |5 lf- SE 28 w*binday — (MT) wrie — mp) |, (S14) 
i 1,9,k,a,8,7 


where e€ is a positive infinitesimal quantity introduced so as to avoid singularities which occur when w is equal to one 
of the eigenvalues of m. To carry out the average over the ensemble of random matrices we need to evaluate 


7 N N aß aß _ Ba 


x det [(yw* — m”)(1yw-m)|] . (S15) 


We next perform a Hubbard-Stratonovich transformation [55] in order to linearise the terms involving meP . This 
yields 


dz ady, 
det [(1yw* — m”) (yw — “= /T oo 573 £) exp - > ye ys + ezz | 


x exp |i 5 gen (mie —w *Sijbap) vf : | exp a> yr * (mee — wõijõap) A (S16) 


ijaß ijaß 


Average over the random matrix and introduction of order parameters 


Carrying out the average over the random matrix elements in Eq. (S16) then gives [recalling the definition of the 
potential ®(w) in Eq. (S8)] 


d2z Pay; 
exp [- w= f Lo E) exp Dur + cite 


x exp Dee + da) + zF yF" (w + da) 


ia 


x exp |—— Y eY HeY +Tapleh yh + eey (2h yt + put). (S17) 


ie [45, 52] we disregard Paar containing N! D ue, Ny (ae) NS yee, See 
-15 2e*yo* and NT! 3°, z*y*. These do not contribute to leading order in N. 
ae now introduce the cider ee 
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Wa = — ayr w* = __ yo” 2 (S18) 


which we impose in the integral Eq. (S17) using Dirac delta functions in their complex exponential representation. 
For example 


1 
A ~: 22 22) x [dies io — Darn] ; 


lwa — —— X 2P*yF) x fe Wa EXP fiaz (NaWa — 5 z? ye) + idal Nawi — 5 ZATE (S19) 


i i 


where we note that wa is a complex quantity. We can thus rewrite Eq. (S17) as 
TR w= fDi D[---]exp(N(W+049)], ($20) 


where where D |---] denotes integration over all of the order parameters and their conjugate (‘hatted’) variables, and 
where 


V= itu + iv tid ya(daws + Wiwa), 
a 


1 
O = —eu — v — o°uv + >D Ya |[—tWa(w* + da) — iw (w + da)] ao 5 Yael og (Wows + wgw), 


a aß 


d?z “dy a 
Q= Dain CE L) exp{ ilz 2% + by y™ + Way z “taney (S21) 


We recall that a can take the values a = u (prey) and a = v (predators), and that yu and 7 are the fraction of prey 
and predator species respectively. 

We note that the integrals over y® and z% are uncoupled for different values of į as a result of introducing the order 
parameters in Eq. (S18). Carrying out the integrals over the variables y“ and z“ in the expression for 2 one obtains 


=X ya naih — ûô). (S22) 


a 


Saddle-point integration and evalulation of order parameters 


We now take the limit N — oo, and carry out the integral in Eq. (S20) in the saddle-point approximation. To do 
this, we extremise the expression Y + © + Q. Extremising with respect to the conjugate variables û, ô, wo, and ù$, 
we find 


$ —YaÛ . —Yaù : Yo y Yo 
1 ` — mi Vz 5 —"*— apn, Wa = Se oe an L a iyaw% = — 22 T (S23) 
One makes the following observation 


itu + idv +iŅ yalÛawi + oh we) = 2. (S24) 


Further, we find that 


x _ Üh 
WaWa = (aÑ = tid)?” (S25) 


which means that we can now express WqwW%, in terms of waw%, and tid. We also find that 


Yo tv 
-w=> $26 
£4 aia — 86) (tpt — GB)’ vey 
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which means that we can in turn express @ in terms of waw%* and uv. Defining r = uv, we can write Gd and awh, 


in terms of ww? and r. This means that we can replace w,w% — ûô in the expression for Q. 
So, letting ga = (waw* — G6)~', we obtain an implicit equation for ĝa in terms of r and waw% : 


We now extremise Ų + © + Q in Eq. (S20) with respect to u,v, Wa and w. We find 
ti=e+o*v, i6=1+07u, 


ide = ilw + da) + 0° DS Tagypwh, it = i(w* +da) +07 5 Top yews: 
B B 


Combining these with Eqs. (S23), one sees that 


v = (e+o7v) a Yoda; 


Q 


u = (1 + 0°u) 5 Yoda: 


Wa = — | ilw +da) +0’ 5 Tag yswp | ĝa, 
B 
wr = — |i(w* +da) +o? S “Tas yews ĝa. 
B 


The first equation in this set has two solutions for €e > 0*. 


First solution: 
One solution is v = 0, implying r = 0, and ĝa = —waw%. The last two of the above equations yield 


1=i(w+da)w, +0? S Tas yews, 
B 
1 = i(w* + da)We + 0? 5 TosyeWaws- 
B 


(S27) 


(S28) 


(S29) 


(S30) 


From this, one can solve for the resolvent G(w) = X „a Yaw%, which we see is an analytic function of w. This means 


that the eigenvalue density vanishes in regions of the complex plane for which r = 0 is the only valid solution. 


Second solution: 
The other solution to the first of Eq. (S29) is 


R 1 
Dafa T F2 


(S31) 


This can be solved simultaneously with the following equations in order to find r, wa and w% as functions of w = z +iy 


Wa = — |ilw +da) +0°S“Tosyews devs 
B 
we = — ji(w* + da) +0°S°Tapyewe da: 


B 


(S32) 
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Significance of the order parameters 


It is important at this point to note that the order parameters defined in Eq. (S18) have some broader significance. 
The quantities wa and wx, when evaluated at the saddle point, are related to the resolvent: 


OP(w,w . A 
Glu) = POH) LD jaut, 
ý aP (w, w* . 
G (w) = (ens = > YaWa- (S33) 


Imagine momentarily that instead of evaluating —N®(a,y) = Indet [(1yw* — m™)(1)w — m)], we had considered 
aß 
tj 
wily in the definition of the resolvent Eq. (S1) with a diagonal matrix (w)°” 


a 
at Eq. (S33), we would have been able to define 


the ‘potential’ of the matrix mi: — wadagdi; [i-e. if we had redone the calculation up until this point but had replaced 


= Wedee6ij]. Then, instead of arriving 


«yy — 1 O®(wa, wah) x 
Ga({Wa, We }) aa Ya dwa = Was 
1 Pwa wat) 


Ya Owes, 


Gi (twas wat) = 


= iwa. (S34) 


The result of this line of reasoning is that, at the saddle point, we deduce iwa = (iw%)*. Such a deduction is still 
valid once we set wa = w. This observation will prove important for ruling out invalid solutions when we solve the 
saddle-point equations. 

Although the fact that at the saddle-point iwa = (¢w*)* may at first seem troubling, one should note the following. 
Let Wa = da + iSo and wx = qa — İSa. The (real) variables over which we integrate in Eq. (S20), in order to impose 
the delta function in Eq. (S19), are qa and Sa. Because the argument of the exponential in Eq. (S20) is complex, we 
have to use analytical continuation to deform the integrals, which are initially over the real axis, to contours in the 
complex plane. These contours include the saddle point. As such, the ostensibly ‘real’ variables qa and sq take on 
complex values at the saddle point. Hence, wa and w% are no longer necessarily complex conjugates when evaluated 
at the saddle point. 


Calculation of boundary of bulk spectrum 


Inside the support of the eigenvalue spectrum (where the eigenvalue density is non-zero), we have r > 0 and the 
condition Eq. (S31). Assuming no discontinuities in the value of r as a function of w, the boundary of the support of 
the eigenvalue spectrum is defined by points w in the complex plane where r approaches zero, but where Eq. (S31) 
is still satisfied. This yields the following set of equations which one can solve for the curve defining the boundary 
Wy(We) (where w = wg + iwy): 


Wa = — | i(we + iwy + da) + o? S “Tas yews ĝa, 
B 

we = — Ji(we — iwy + da) + g? S “Tas yews ĝa- (S35) 
B 


That is, values of w which lie on the boundary of the support of the eigenvalue spectrum must satisfy the above 
simultaneous equations. 

Since iw, = (tw*)* at the saddle point, the second relation in Eq. (535) dictates that ĝa is real and positive. In 
solving Eqs. (S35), we therefore discard solutions which do not satisfy ĝe € R and ĝa > 0. Alternatively, letting 
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iw = (iwa)* = Xa, we obtain 
` Yo Xal : ’ 
a2 
Q 


1 = — (wr + iwy 4 da)Xa +0? X` TagY8XaXp: (S36) 
B 


S3. OUTLIER EIGENVALUES 
A. Setup and general relations 


We now wish to find the outliers that arise when we perturb the random matrix m by introducing a finite mean. 
That is, we are now interested in the matrix m’ = m+ x H, where the NxN matrix u accounts for the non-zero 
mean. Specifically m’ has elements 


mP = meP q PR (S1) 
in particular 


(u)$f = pap, (82) 


independently of i and j. The elements of m have zero mean as before. So, the elements within each of the 4 
blocks [(a, 8) = (u, u), (u,v), (v, u), (v, v)] of the matrix m’ have the same mean, pag/N. Following the conventions 
of statistical physics, the mean is chosen to be proportional to 1/N for technical convenience. In the main paper, 
the mean is independent of N. We explain the mapping between the two setups in Section $4, where we also discuss 
general values of connectance 0 < C < 1. 

The problem amounts to finding the eigenvalues of the matrix m’ — i.e. m with a rank-2 perturbation u. We follow 
a procedure similar to Ref. [47]. The bulk of the eigenvalue spectrum is unchanged by the low-rank perturbation. 
Therefore we look for additional eigenvalues that lie outside of the bulk of eigenvalues of m. Such an eigenvalue A 
must satisfy 


1 
det [at —m-— A =0. (S3) 


Any outlier A, is by definition, not an eigenvalue of the unperturbed matrix m. Therefore the matrix 1yA — m is 
invertible. Hence, 


det fix- ly — my =0. (S4) 


B. Simplification to an effective 2-species problem 


Now we define the resolvent matrix G = (Aly — m)7?. 


Case Nu = Ny: 

We suppose for now that there are the same numbers of prey species as predator species, i.e. Na, = N,. We will relax 
this assumption once we have outlined the proof for this simpler case. Suppose we rearrange the rows and columns 
of the matrices so that the first species is a predator species, the second a prey species, the third again a predator 
species, and so on. The matrix p then consists of N, x N, block of size 2 x 2 of the form p?) = ( Huu Bean I Then, 


Hvu Luv 
the matrix u can be written in the form 


=u, (55) 


where u and v are vectors of length N,,. Each entry of v is given by the 2 x 2 matrix p®, 


ges (un, a -) (S6) 
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The vector u is of the form 
u = (lz, 12, l2,- )7. (S7) 
We next use Sylvester’s determinant identity 
det [Im + A B] = det [1 + BA], (S8) 


valid for combinations of m x k matrices A and k x m matrices B. We note that the matrices on the left-hand side 
of Eq. (S8) are of size m x m, and those on the right-hand-side of size k x k. 
Using this identity, we find that det(Ilz — v7 (vA — m)~!u) = 0.Thus we obtain 


1 E 
eyo) ae (ij) _ 
det | 12 — u a. =0, (S9) 


= ((lyA—m)-1)2? . 
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where (GED) að 


General case (N, not necessarily equal to N»): 
We now wish to show that Eq. (S9) also applies in the case Nu # N,. This is done using the following identity for 
block matrices: If A, B, C and D are matrices of sizes m x m, m x k, k x m and k x k respectively, then we have 


AB A 0 
det E J = det(A) det(D) = E J : 

A B] [A B] _ [AA’+ BC! AB’ +BD' 
C D| |œ” D CA + DC’ CB'+ DD'|: 


(S10) 


First we note that the determinant det |1 [il N— + Gu] in Eq. (S4) is unchanged by adding extra rows and columns 
to the matrices in the following way: We add rows and columns yp to produce a square matrix fy. of dimension 
N* = N + |N, — N,| which can be decomposed as in Eq. (S5). We then augment G to obtain 


G 0 
ex- [2 4) m 
Then, using the theorems in Eq. (S10), one can show 
det | 1 ld = det |1 la (S12) 
D ee | ee ee 


Hence, we can apply Sylvester’s determinant identity as before and arrive at Eq. (S9). 


C. Proof that the resolvent matrix is diagonal 


We have shown that Eq. (S9) holds in all cases. It remains now to find the averaged Green’s functions 
+ Dij GI) (A). For this, we follow [51] and use the so-called ‘locator expansion’ to show that the resolvent ma- 


trix is diagonal in the limit N — oo, i.e. that Ge? x dijðap. Letting Gi = (f-") 608615 — mee (where in our 
case (f—~+)% = A), it follows that m = G-1(G — Pf- r iom which we obtain a Dyson equation G = f + Gmf. This 
allows us to construct a series for the resolvent 


ô 
= fE HIEP SE + MERNE + famey RM mP f +... , (S13) 


where we have ae sums over repeated indices for brevity. We find that we have to evaluate terms such as 


— ay yB PER ; ; ; 
Zorg m 5 aTe mis z5 ka fima Myj- This involves averages over products of Gaussian variables. Such terms can be 
simplified using ’ standard Gaussian combinatorics 


m21! m2222 co. monn — — 5 II meee, (S14) 


it izi inih 
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which is valid for even n and is zero otherwise, and where }~ |[ stands for the sum of the products of all possible 


Qaba R 


combinations of pairs m; cint M; Since we have 
Qaba Abbo g? Toata 
Miis Migi = Py iei 914,74, aaar 98.85 + Sig it, Oi, 1 i Seva Br 9Bacs (S15) 


one can see that (since the only free indices are i, j, a and 6 in Eq. (S13)) that Gee xX 6i;5a8- We therefore find that 
2G GRP (A) © bap 2 Gy? (S16) 


for large N. 


D. Final expression for outliers 


Given that the resolvent matrix is approximately diagonal for large N, we only need to calculate the sum of diagonal 
elements 4 X; G&* in Eq. (S9). This can be found from the quantities involved in the bulk calculation as follows. 


Consider the potential ®({w.,w%}) associated with the generalised resolvent G({wa}) = + oi, (w — m)~!, where 
(w)? = = Waba86i;- We know that at the saddle point 


— pt a 1 llw) 


= 1W 
Xa a Ya Owe, 


(S17) 


The potential is defined by ®({wa,w%}) = —#Indet [(w — m)t (w — m)]. Now, differentiating with respect to wa and 
using Jacobi’s formula for the derivative of a determinant, we have 


Pwa, wx }) 1 
Bu = wit c m)~ 


: 2e- = | = 5 WD G2" ias): (S18) 


Now setting wa = w, we see from Eq. (S9) that the outliers À obey 


det [1> — p” yx0)] =0, 


- Yas 0 _ |Xu 0 
= |e 2], xar= pe 2]. (S19) 
This can be written as 
| am | = | =0 ($20) 
YuHuu Xu A) Yu Huv Xo A) Yu Yv Huv Hvu . 


The quantities Xa(A) satisfy [c.f. Eq. (S36)] 


—1=-(A + du) Xu + oT aux? + oT uv teXoXur 
-l= —(A Bg dy) Xv + OD uw YuXuXe + gl Wie oe (S21) 


Note that not all solutions of this set of equations are valid. At the saddle point, the order parameter u = x yee A 
must satisfy u > 0, and we have u = (1+07u) >>, YalXo|?. This means that we must have 


Enka? < = (822) 


in order for the set (A, {xa(A)}) to be a valid solution to the combination of Equations (S20) and (S21). 
Eqs. (S20) and (S21) can be solved numerically for A, subject to the constraint in (S22). Note that when the 
left-hand side is equal to the right-hand side in the inequality (S22) the point A in question is inside the bulk region. 
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S4. SCALING WITH N AND CONNECTANCE C 


There are two key differences between the random matrices used in Sections S2 and $3 [c.f. Eq. (S12)] and those in 
the main text [see Methods section]. These are: (a) The model parameters are rescaled by factors of N; (b) Not all 
species interact with all others in the main text (i.e. the connectance can take values smaller than one, 0 < C < 1). 
We now describe how to relate the results of the previous sections to the setup in the main text. 


A. Rescaling with N 


In this Supplement, we used the statistics given in Eqs. (S12) and (S1) for the entries of the random matrix m. 
These are 


aß _ Haß 
my =N 
aß o? 
(mi; )? = (Hap) = N’ 
TBT T B 1107 
(mg — pap (ngs — pap) = OO (S1) 


The statistics used in Eq. (4) and (5) of the the main paper (and also given in Eq. (S2) of this Supplement) are 
different. They can be obtained from the above by making the replacements pag > Nuag and o? + No?. The 
results in the main paper have been obtained from those in Sections $2, S3 and S5 by making this replacement. 

The theoretical predictions for stability or instability are formally derived in the limit N — oo, but in practice 
they are accurate for finite values of N, provided N is not too small. Sizes of N ~ 50 are often sufficient to obtain 
satisfactory agreement (see Fig. S1 below for an example). The spectra in Fig. 2 of the main paper are obtained 
numerically from matrices of size 750 x 750, and can be seen to agree very well with the results of the mathematical 
theory. 


B. Taking into account connectance C 


A connectance C < 1 between species can be accounted for by making a similar replacement of o? and pag to the 
one described above. We now describe how to see this. 

In order to find the potential ®(w) (which contains all the information about the eigenvalue spectrum and could 
in principle also be used to find the outliers [52]) defined in Eq. (S8), one needs to carry out an average over the 
ensemble of random matrices, see e.g. the step from Eq. (S16) to Eq. (S17). We now wish to take into account that 
any one pair of species only interacts with probability C. We then obtain [c.f. Eq. (S15)] 


I II (amf) P({m$f}) det [(Lnw* — m?) (yw — m)] 


ijaß 


= f TY (ame?) det [Cwe - m?)(a,yee = m)] 


ijaß 


y [f — Hap)? = Palmi? — Has) mie — vo) \ (82) 
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Upon evaluating the integral, one obtains 
d? zl d2ye 

exp [—N®(w)] = Ju (=) exp - a + ate 

x exp -: NO eye (w* + da) + zry (w + da) 


x II {u=e)+ Cex 


i 
ples (2 "yf + 2p vs") 


o2 
— noa HEYS + Dap (2htyh + eyf ei ye + zBye*)| a. (S3) 


If we set C = 1, we recover the term we would normally obtain with all species connected 
dz% d? y% 
exp [—N®(w)] = JU (E) exp | — Daas + Eze Zo 
ia ia 


x exp | -iX 28 *y%(w* + da) + 2h y?*(w + da) 
ia 


i 
Xx exp È a, Hop (o +z; us") 
i j,a, b 

2 


oO 
-3N De RUG + Pus")? + Papleh ah + aR uy Vey" uF + 2 0] ; (54) 


For individual sets of i, j, a and 6, and assuming N > 1, the expression in the curly brackets in Eq. (S3) can be 
expanded in orders of 1/N to obtain 


i x a ax 
{u-)+ Cex les (a Y7 + zi Yj ) 


2 
oO 
— ay Puy tau V + Daal yy eu Vey Us + 25 w) j 


i 
~ 1+ CHol GURE vs") 
Cor neg a, x2 ax, b a, Bx \7 Bx a Bo ax 
— Sar ee Wg tau y Hagley + 2iyy G ue + 25 92") 


a 
= exp gCo (yg + z; us") 
Co? ax, B a, Bx\2 ax, B a, Bx Bea Bak 
-VN (zi u Heu y + Taglar yy tziu G y tzu) (S5) 


Comparing Eq. (S3) [with Eq. (S5) substituted in] to Eq. (54), we see that the consequence of introducing a connectance 
0< C <1is to scale the system parameters in the following way 


Hop + Chap, 
o? + Co’. (S6) 


C. Combined effect 


In summary, the scaling that one needs to perform in order to obtain the results in the main text from those in 
Sections $2 and S3 is given by 


Hap > NC uap, 
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o? + NCo*. (S7) 


All simulation results in the main text were performed using the model as defined in Section S1 (which is the model 
defined in the main text). We note that Eqs. (S36) and (S20-S22) become respectively 


yt Xal — : = 0, 


—(Wy + iwy + da + g Da)Xa + c> Tap yexaXs +1=0, (S8) 
B 
and 
1 1 
YuN C huu = Yo NC how om Nl (NC) YuYu Huv Hou = 0, 
Xu(w) Xv(w) 


—(w +dy+ g Du)Xu + Tuu Xi + CPuyYuXuXu =-l, 
—(w + dy + Dy) Xo + VuvYuXuxe + Poyxy = —1, 


SW relxale)l? < =, (89) 


where c = NCo?. 


S5. STABILITY CRITERIA 


We class the instabilities in our system into two types: 


(1) The bulk of the eigenvalue distribution may cross the imaginary axis. This is the type of instability May 
identified [7]. 


(2) One of the outlier eigenvalues may cross the imaginary axis, with the bulk staying in the negative half-plane. 
Instabilities of this type are discussed for non-spatial systems for example in [13]. 


Below, we provide analytical criteria for each of these instabilities. Crucially, we note that an instability of the second 
kind can be introduced to an otherwise stable system by the inclusion of diffusion. That is, for a finite band of 
wavenumbers, the spatial system is unstable where the non-spatial system would be stable. 

First, we provide criteria for the stability of the non-spatial system, and then we characterise the spatial instability. 
In the following the model parameters o° and Hag are scaled in the same way as in the main text, and we allow for 
general connectance 0 < C < 1. 


A. Instabilities of the non-spatial system 
Instability due to bulk eigenvalues, and relation to May’s bound 


Mathematical condition for instability: 

The transition occurs when the boundary of the bulk eigenvalue spectrum first hits the imaginary axis. We generally 
observe that the bulk of the eigenvalue spectrum is a convex set. Since the matrix we are dealing with is real, we 
also know that the spectrum is symmetric with respect to the real axis. As a consequence we find that point of first 
contact is at w = 0. We note that, from the correspondence with the resolvent shown in the previous section, this 
means that the variables Xa(0) must be real. In order for the system to be stable, A = 0 must necessarily be outside 
the bulk (otherwise there are positive real eigenvalues). Hence, a necessary condition for stability is 


1 
Do (S1) 


We have used Eq. (S22) with the rescaling pag > C'Hap, and o? —> Co?, and we recall that the complexity parameter 
is given by c = NCo?. The quantities Xa satisfy 


—1 = —daXa + cX Tap exaXxe- (S2) 
B 
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The equations above cannot be solved analytically in general, but a numerical solution can be obtained. 


Special cases: 

We can make the connection with May’s original criterion [7] in several special cases. First, consider the case where 
Yy = 0 — that is, we do not distinguish between predator and prey species, and consider a general population of 
N = N, species (as May did). In this case, Eqs. (S1) and (S2) yield 


Ve(1+Tuu) < du. (S3) 


We recover May’s criterion for Tuu = 0 and du = 1. Secondly, the criterion for stability also simplifies greatly when 
the prey and predator populations have uncorrelated internal interactions so that Tą =I, = 0, the same number of 
species (i.€., Yu = Yv = 4), and equal death rates du = d, = d. In this case, Eqs. (S1) and (S2) give 


1 
v2c(1 + grw) <d. (S4) 
Similarly, if Tù = 0 and Tu = Te =T, we obtain 


Vli Ir) <d. ($5) 


A further special case is a situation in which all interactions are uncorrelated (Tag = 0 for all a and 8). Here we find 


ey 225, (S6) 


From these examples, we get a general idea of how each of the model parameters affect the onset of instability. The 
intra-species interaction Ia = —da is required to be negative in order to have stability. The more negative Ia is, the 
more stable the equilibrium is. The equilibrium can also be stabilised by asymmetry in the interaction coefficients. 
That is, the lower the parameters —1 < Iag < 1 are, the more stable the equilibrium. Finally and crucially, as per 
May’s conclusion, a lower complexity c = NC? contributes to making the equilibrium more stable in all the special 
cases outlined above. 


Instability due to outlier eigenvalues 


Mathematical condition for instability: 

If we suppose that the bulk eigenvalues remain to the left of the imaginary axis, an instability may still occur if one or 
more of the outlier eigenvalues take on a positive real part. This may occur in two possible ways: (I) A single outlier 
with no imaginary part may cross at A = 0; or (II) A pair of complex conjugate outliers may cross the imaginary 
axis. The point at which an instability of the first kind occurs can be obtained by substituting A = 0 in Eqs. (S20) 
and (S21). Re-scaling o? and fag to match the conventions in the main paper, we obtain 


1 
2 = (NC)? Yu Yv bur bow =0, 


v 


1 
[NOus = + [Cts = 


uU 


Vuyux2 + CPunYuXuXu = duXu +1=0, 


CDuvYuXuxv + Towi = dyXv +1=0. (S7) 


The model parameters at which an instability of type (II) above occurs are found as follows: We evaluate the outlier 
eigenvalues by solving Eqs. (S20) and (S21) subject to Eq. (S22), and we carry out the re-scaling pag > NC Hap, 
a? — NCo?. We do this for varying model parameters until one of the outliers gains a positive real part. This way, 
we ensure that we spot any complex conjugate pairs crossing the imaginary axis. The sloped lines in Fig. 3 in the 
main text, for example, were generated in this way. 


Example: 

Fig. 2 in the main text shows an example in which a single outlier crosses the imaginary axis. In Fig. S1 we show 
a case in which a complex-conjugate pair of outliers has crossed the axis. The model parameters are the same as in 
Fig. 2 of the main text, with the exception of the uag. This shows that altering the means of the matrix elements 
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can turn an instability of type (I) above into one of type (II). We note that the bulk of the eigenvalue spectrum is the 
same in Fig. 2a and in Fig. $1. This confirms that the bulk spectrum is not affected by changes in the means Hag. 
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FIG. S1. Example eigenvalue spectrum showing a pair of eigenvalues crossing the imaginary axis. Eigenvalues obtained 


numerically from computer-generated random matrices are shown as blue crosses. Red lines are the theoretical predictions for 
the boundary for the bulk region, and red hollow circles the predictions for the locations of any outliers. The model parameters 
are the same as in Fig. 2 (a) (with q = 0) of the main text, with the exception of the means pag which are now given by 
CN huu = 4, CN buv = —4, CN hou = 4, CN bow = —2, where the scaling of the system parameters is the same as in the main 
text. In the left-hand panel we use N = 750, on the right N = 45, keeping the products CN uag constant. 


Special cases: 
We can make analytical progress in the same set of special cases as in the previous section. First, consider the case 
Yu = 0 (we could similarly use yẹ = 0). Then, the outlier eigenvalue [from Eq. (S21)] is real and given by 


Tyc 
A= NC loo S8 
How + Na, (58) 
For C = 1 this reduces further to results mentioned in [47, 52]. A necessary condition for stability is then 
Tye 
NC hov + —— <0. 59 
bo + ROm < (59) 


In the case where the prey and predator populations have uncorrelated internal interactions so that Tą =T, = 0, 
the same numbers of species yy = Ww = 3, and the same death rates d, = dy = d and Huu = —Mvyv, the outlier 
eigenvalues are given by 


N 2 2 UVP UU 2w 
y d4 WO an + Hutton) + Pwe (S10) 


2NC V Hu + Bulow 


If p2.,, + Mutou < 0, there is no real solution for \ and there is guaranteed to be no instability, so long as d > 0. We 
note that complex conjugate pairs of outliers are made possible only by having different trophic levels. 
Similarly, in the case Tus = 0 and Tu = Te =, the eigenvalues are given by 


N 2 UV PUUW 2r 
y q4 WO Hiu + Huvhvu) + Pc (S11) 


2NCV/UZu + Mav bow 


In both special cases, the more negative the quantity HuvHvu is, the more stable the equilibrium. We note that 
— HuvHvu is a measure for the ‘strength’ of predation (taking into account both the effects on prey and predators). 
So, predation is a stabilising influence. Further, asymmetry in the interaction coefficients has a stabilising effect. 

A further special case is Tag = 0 for all œ and 8. Here, we obtain 


[NC Yubbuu = du] [NC Vv buy _ da] = (NO)? Yu Yu Hav ton > 0. (S12) 


Once again, the more negative the quantity HuvHvu is, the more stable the equilibrium. 
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B. Dispersal-induced instability 


So far we have considered only a non-spatial ecosystem in this Supplement, i.e. we have set q = 0 in Eq. (S4), or 
equivalently D, = D, = 0. We now address instabilities induced by dispersal. The additional terms arising from a 
non-zero dispersal term can be absorbed into a re-definition of da, 


da > da oe Ds (S13) 


This can be seen directly from Eq. (S6). With the replacement in (S13), we can then use the stability criteria derived 
in Sec. S5 to obtain criteria for dispersal-induced instability. 


The instability due to bulk eigenvectors is unaffected by dispersal 


With the inclusion of dispersal Eqs. (S1) and (S2) become 
1 
D eee 
z c 


—1 = —(da +° Da)Xa + cX Top yeXaXs- (S14) 
B 


This is the condition for the stability of the Fourier mode with wavenumber q. 

Although it is difficult to prove analytically, we find in general that the inclusion of dispersal serves only to shift 
the bulk of the eigenvalues in the negative real direction, never in the positive. This is seen to be the case in Fig. 
2 in the main text. This can also be seen to be true in the special cases we have discussed already — increasing the 
magnitude of da (which is effective what one does by introducing non-zero q) in Eqs. (S3)-(S6) serves only to expand 
the range of parameters for which the equilibrium is stable. 


Turing-type instability due to outlier eigenvalues 


The Turing instability is a ‘stationary’ instability [56] which exhibits no oscillations. This means that the eigenvalue 
crossing the imaginary axis has no imaginary part (unlike the example given in Fig. S1). The point in parameter 
space at which an instability of this type occurs in the non-spatial case (q = 0) is given by the solution of Eq. (S7). 
Using this, and the replacement da > da + q?D, to account for q Æ 0, we arrive at the condition for stability with 
respect to disturbances of non-zero wavenumber q: 


1 1 
P(q) = [Catan > al [Coto = > ~~ (NC) Wire tuclten > 0, 


Duyux2 + PuvYXvXu — (du + Dud?) Xu +1=0, 
CPuvYuXuXv T Tuvi = (dy + Dad" Xv +1=0. (S15) 


In order for the ecosystem to show a Turing instability this condition must be satisfied for q = 0 (i.e., the system is 
stable with respect to perturbations with q = 0), but violated for some non-zero q (there are unstable modes with 
q #0). That this can indeed be possible is exemplified by the special case where Iag = 0 for all a and 8. In this 
case, the stability criterion for the Fourier mode q becomes 


Polg’) = [NC yu tun — (du + Dig’) [NC yb — (dy + D,@)| = (NO) Poe) Papa > 0. (S16) 
The expression for Po(q) is quadratic in q?. It has a minimum at 


N uuu ~ u)Dy + (N v vy — dy) Du 
ae Cuh ta a Cw Du (817) 


If Po(0) > 0 and Po(q2\;,) < 0, then there is a Turing instability. This occurs when 


Dy 1 
—— Yu Yu 7 
YuHuu 


2 
D. (i — Huvhvu + VEFi) , (S18) 
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where NC 7 ut, = NC Yubtuu — du and NC wp), = NCh — dy. This is very reminiscent of the criterion for 
Turing instability in more conventional reaction-diffusion systems [27]. 

In the more general case, we need to find the minimum possible value of the expression on the left-hand side of the 
first relation in Eq. (S15). If the minimum value of P(q?) is negative (i.e. the condition in the first line of Eq. (S15) 
is violated), but the value of P(0) is positive, then we say that there is a Turing (dispersal-induced) instability for 
this parameter set. 

Again, we first find the minimum value of P(q?) for a given set of model parameters. To do this, we differentiate 
P(q’) with respect to q?, noting that the three relations in Eq. (S15) are coupled. This leads to the following set of 
simultaneous equations to find P(q2in) 


P' (Qin) = NCYutuuxy (NC WwXvtivv — 1] + NC yt xy [NCyuXubun — 1] 
— (NC) Yu Yottuvtivn XuXy + Xoxa] = 0, 
O = WD ua Vu Xu Xa + Pave (XoXua + XoXu) — (du + Duqin) Xa — DuXu» 
0 = Duna (XoXa + XXu) + 2M vw WwXvxy — (dy + DuGhin)Xy — Dox: 
O = ye, + PuvWwXvXu — (du + Duta) te + 1, 
0 = Dav tuxuxe + Vox — (dv + Dodhin)Xo + 1, (S19) 


where y/, = Xa, One solves the above equations simultaneously to find yu, X4, Xv; Xj, and qŽin: One then evaluates 


PG) = [NC —4uXukHuu = 1] [NC wXvhvv = 1] T (NC) YuYohuvhvuXuXv: (S20) 


So, one can find the minimum value of P(q?) for a particular parameter set and thus deduce whether or not there 
is dispersal-induced instability (i.e. whether the condition in Eq. (S15) is satisfied). For the purpose of constructing 
stability diagrams, one can instead set P(q?,;,,) = 0, equivalent to assuming that the outlier eigenvalue is just about to 
cross the imaginary axis (A = 0). One then solves for the critical model parameters instead which give P(q2in) = 0. 
This can be used to produce solid sloped lines in Figs. 3 (b) and 4 in the main text . 

Although the above set of equations appears cumbersome, we note that they are at most linear in x/, and x/,, so 
these can be solved for relatively easily. In practice, Eqs. (S19) reduce to three simultaneous equations for qin; Xu 
and Xv. 


C. Dispersal cannot make an unstable equilibrium stable 


In our model there are no instances in which an unstable equilibrium of the non-spatial ecosystem becomes stable 
when dispersal is introduced. This is because a stable equilibrium in the spatial model must be stable with respect 
to disturbances for values of q, including q = 0. This means that such an equilibrium is also stable in the non-spatial 
system. 


S6. PARAMETERS USED FOR FIGURES IN THE MAIN TEXT 


Figure 2: 
N = 750, C = 0.5, yu = 2/3, w = 1/3, c = CNo? = 0.25, du = 1.0, dy = 16,71, = 0.5, Tuv = —0.9, rẹ = —0.5, 
CN buu = 2, CN bws = —2, CN tvu = 2, CN ktos = —1, Dy = 1, Dy = 5. 


Figure 3: 
T, = 0.5, Ty = —0.9, Puy = —0.5, Dy = 1, Dy = 5, dy = dy = 1, CN puu = 2, CN pos = —1, Yu = 2/3, w = 1/3, 
Huv = —Hvu. 


Figure 4: 
Pa = 1.0, P, = —1.0, Pus 0, Du 1, dy dy 1, CN huu = CN hou =2, CN uw = —2, CN bow =-l, Yu = 2/3, 
äp = 1/3. 
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Figure 5: 
N = 100, T, = 0.5, Py = —0.5, Py, = 0, Dy = 0.2, Dy = 20, VC No = 0.1, dy = dy = 0.1, CN huu = 0.6, CN hvu = 1, 
(b) Dy 


CN tuw = —1, CN hoy = —0.4, a = 0.5, Yu = 0.5, w = 0.5. In (a), Du, Dy = 0 and in (b) 


= 0.2, D, = 20. 


Eugene P. Odum and Gary W Barrett, Fundamentals of ecology (Saunders, Philadelphia, 1953). 

Charles Sutherland Elton et al., Ecology of invasions by animals and plants (Wiley, London, 1958). 

Robert MacArthur, “Fluctuations of animal populations and a measure of community stability,” Ecology 36, 533-536 
(1955). 

Robert T Paine, “Food web complexity and species diversity,” The American Naturalist 100, 65-75 (1966). 

Pietro Landi, Henintsoa O Minoarivelo, Ake Brannstrém, Cang Hui, and Ulf Dieckmann, “Complexity and stability of 
ecological networks: a review of the theory,” Population Ecology 60, 319-345 (2018). 

Kevin Shear McCann, “The diversity—stability debate,” Nature 405, 228-233 (2000). 

Robert M May, “Will a large complex system be stable?” Nature 238, 413-414 (1972). 

Robert M. May, “Stability in multispecies community models,” Mathematical Biosciences 12, 59 — 79 (1971). 

Toshiyuki Namba, “Multi-faceted approaches toward unravelling complex ecological networks,” Population Ecology 57, 
3-19 (2015). 

James Justus, “A case study in concept determination: Ecological diversity,” in Philosophy of Ecology, Handbook of 
the Philosophy of Science, Vol. 11, edited by Kevin deLaplante, Bryson Brown, and Kent A. Peacock (North-Holland, 
Amsterdam, 2011) pp. 147 — 168. 

Jacopo Grilli, Tim Rogers, and Stefano Allesina, “Modularity and stability in ecological communities,” Nature Commu- 
nications 7, 1-10 (2016). 

Stefano Allesina and Si Tang, “The stability-complexity relationship at age 40: a random matrix perspective,” Population 
Ecology 57, 63-75 (2015). 

Stefano Allesina and Si Tang, “Stability criteria for complex ecosystems,” Nature 483, 205-208 (2012). 

Lewi Stone, “The feasibility and stability of large complex biological networks: a random matrix approach,” Scientific 
Reports 8, 1-12 (2018). 

Theo Gibbs, Jacopo Grilli, Tim Rogers, and Stefano Allesina, “Effect of population abundances on the stability of large 
random ecosystems,” Phys. Rev. E 98, 022410 (2018). 

Eric L Berlow, Anje-Margiet Neutel, Joel E Cohen, Peter C De Ruiter, BO Ebenman, Mark Emmerson, Jeremy W 
Fox, Vincent AA Jansen, J Iwan Jones, Giorgos D Kokkoris, et al., “Interaction strengths in food webs: issues and 
opportunities,” Journal of Animal Ecology 73, 585-598 (2004). 

Thilo Gross, Lars Rudolf, Simon A Levin, and Ulf Dieckmann, “Generalized models reveal stabilizing factors in food 
webs,” Science 325, 747-750 (2009). 

Peter Chesson and Nancy Huntly, “The roles of harsh and fluctuating conditions in the dynamics of ecological communities,” 
The American Naturalist 150, 519-553 (1997). 

Samuel J McNaughton, “Diversity and stability of ecological communities: a comment on the role of empiricism in ecology,” 
The American Naturalist 111, 515-525 (1977). 

Elisa Thébault and Michel Loreau, “Trophic interactions and the relationship between species diversity and ecosystem 
stability.” The American Naturalist 166, E95-E114 (2005). 

David Tilman, “The ecological consequences of changes in biodiversity: a search for general principles,” Ecology 80, 
1455-1474 (1999). 

Alan Mathison Turing, “The chemical basis of morphogenesis,” Bulletin of Mathematical Biology 52, 153-197 (1990). 
Ilkka Hanski, “Metapopulation dynamics,” Nature 396, 41—49 (1998). 

Michael P Hassell, Hugh N Comins, and Robert M May, “Spatial structure and chaos in insect population dynamics,” 
Nature 353, 255-258 (1991). 

Michael P Hassell, Hugh N Comins, and Robert M May, “Species coexistence and self-organizing spatial dynamics,” 
Nature 370, 290-292 (1994). 

Simon A Levin and Lee A Segel, “Hypothesis for origin of planktonic patchiness,” Nature 259, 659 (1976). 

Michael Cross and Pierre Hohenberg, “Pattern formation outside of equilibrium,” Rev. Mod. Phys. 65, 851-1112 (1993). 
Michael G Neubert, Mark Kot, and Mark A Lewis, “Dispersal and pattern formation in a discrete-time predator-prey 
model,” Theoretical Population Biology 48, 7-43 (1995). 

Max Rietkerk and Johan Van de Koppel, “Regular pattern formation in real ecosystems,” Trends in Ecology & Evolution 
23, 169-175 (2008). 

Reinier HilleRisLambers, Max Rietkerk, Frank van den Bosch, Herbert HT Prins, and Hans de Kroon, “Vegetation pattern 
formation in semi-arid grazing systems,” Ecology 82, 50-61 (2001). 

Simon Levin, “Dispersion and population interactions,” The American Naturalist 108, 207-228 (1974). 

JD Murray, Mathematical biology I: spatial models and biomedical applications (Springer New York, 2001). 

Joseph W. Baron and Tobias Galla, “Stochastic fluctuations and quasipattern formation in reaction-diffusion systems with 
anomalous transport,” Phys. Rev. E 99, 052124 (2019). 

Jacopo Grilli, Matteo Adorisio, Samir Suweis, György Barabás, Jayanth R Banavar, Stefano Allesina, and Amos Maritan, 


35 


36 


37 


38 


39 


40 


41 


42 


43 


44 


45 


46 


47 


48 
49 


50 
5l 


52 


53 
54 


55 
56 


30 


“Feasibility and coexistence of large ecological communities,” Nature Communications 8, 14389 (2017). 

Yoshiki Kuramoto, “Diffusion-induced chaos in reaction systems,” Progress of Theoretical Physics Supplement 64, 346-367 
(1978). 

Mercedes Pascual, “Diffusion-induced chaos in a spatial predator-prey system,” Proceedings of the Royal Society of London. 
Series B: Biological Sciences 251, 1-7 (1993). 

Istvan Lengyel and Irving Epstein, “A chemical approach to designing turing patterns in reaction-diffusion systems.” 
Proceedings of the National Academy of Sciences 89, 3977-3979 (1992). 

Vincent Castets, Etiennette Dulos, Jacques Boissonade, and Patrick De Kepper, “Experimental evidence of a sustained 
standing turing-type nonequilibrium chemical pattern,” Phys. Rev. Lett. 64, 2953-2956 (1990). 

Dominique Gravel, Francois Massol, and Mathew A Leibold, “Stability and complexity in model meta-ecosystems,” Nature 
Communications 7, 12457 (2016). 

Jacob D O’Sullivan, Robert J Knell, and Axel G Rossberg, “Metacommunity-scale biodiversity regulation and the self- 
organised emergence of macroecological patterns,” Ecology Letters 22, 1428-1438 (2019). 

Nicholas J Gotelli, Hideyasu Shimadzu, Maria Dornelas, Brian McGill, Faye Moyes, and Anne E Magurran, “Community- 
level regulation of temporal trends in biodiversity,” Science advances 3, e1700315 (2017). 

Anne E Magurran, Amy E Deacon, Faye Moyes, Hideyasu Shimadzu, Maria Dornelas, Dawn AT Phillip, and Indar W 
Ramnarine, “Divergent biodiversity change within ecosystems,” Proceedings of the National Academy of Sciences 115, 
1843-1847 (2018). 

James H Brown, SK Morgan Ernest, Jennifer M Parody, and John P Haskell, “Regulation of diversity: maintenance of 
species richness in changing environments,” Oecologia 126, 321-332 (2001). 

Jennifer M Parody, Francesca J Cuthbert, and Ethan H Decker, “The effect of 50 years of landscape change on species 
richness and community composition,” Global Ecology and Biogeography 10, 305-313 (2001). 

Fritz Haake, Felix Izrailev, Nils Lehmann, Dirk Saher, and Hans-Jürgen Sommers, “Statistics of complex levels of random 
matrices for decaying systems,” Zeitschrift fiir Physik B Condensed Matter 88, 359-370 (1992). 

Hans Juergen Sommers, Andrea Crisanti, Haim Sompolinsky, and Yaakov Stein, “Spectrum of large random asymmetric 
matrices,” Phys. Rev. Lett. 60, 1895-1898 (1988). 

Sean O’Rourke, David Renfrew, et al., “Low rank perturbations of large elliptic random matrices,” Electronic Journal of 
Probability 19 (2014). 

Endre Süli and David F Mayers, An introduction to numerical analysis (Cambridge university press, 2003). 

Tobias Galla, “Random replicators with asymmetric couplings,” Journal of Physics A: Mathematical and General 39, 3853 
(2006). 

Guy Bunin, “Ecological communities with lotka-volterra dynamics,” Phys. Rev. E 95, 042414 (2017). 

Alan J. Bray and Michael A. Moore, “Evidence for massless modes in the ‘solvable model’ of a spin glass,” Journal of 
Physics ©: Solid State Physics 12, L441—L448 (1979). 

Samuel F Edwards and Raymund C Jones, “The eigenvalue spectrum of a large symmetric random matrix,” Journal of 
Physics A: Mathematical and General 9, 1595 (1976). 

David Sherrington and Scott Kirkpatrick, “Solvable model of a spin-glass,” Physical review letters 35, 1792 (1975). 

Marc Mézard, Giorgio Parisi, and Miguel Virasoro, Spin glass theory and beyond: An Introduction to the Replica Method 
and Its Applications, Vol. 9 (World Scientific Publishing Company, London, 1987). 

John Hubbard, “Calculation of partition functions,” Phys. Rev. Lett. 3, 77-78 (1959). 

A. Yadav and Werner Horsthemke, “Kinetic equations for reaction-subdiffusion systems: Derivation and stability analysis,” 
Phys. Rev. E 74, 066118 (2006). 


